function figure_handle = rt_plotdat(t_interval, t, P, T, S, O, text, av_win, av_thresh)
    % Create the plots
    %
    % (Typical) Usage:
    %
    % Typical usage is that this function is called by rt_driver; 
    % typically, the user will not call this function directly
    
    global siteID

    % Default averaging threshold is 70 samples, window is 1 day
    if nargin < 9
        av_thresh = 70;
    end
    if nargin < 8
        av_win = 1;
    end
    skip_contour = false;
    try
        % Create figure
        figure_handle = figure('units', 'inches',...
                               'color', 'w',...
                               'position', [0 0 8.5 11],...
                               'paperposition', [0.5 0.5 7.5 10]);

        import containers.*
        % Bottom-left corners and heighths and widths of the time series plots
        tserieslp = 0.1; tseriesbp = [0.8 0.65 0.5 0.35];
        tseriesw = 0.80; tseriesh = 0.15;

        % Calculate the starts and ends of the night-indicating patches
        sv = datevec(t_interval(1)); 
        shour = sv(4);
        if (2 < shour) && (shour <= 14) % if start is during the night (locally)
            % First patch begins immediately, subsequent patches begin at hour 2, starting tomorrow
            patchstarts = [t_interval(1) datenum([sv(1:2) sv(3)+1 2 0 0]):1:t_interval(2)];
            % All patches-except perhaps the last, to be dealt w/ in next code block-end at hour 14, 
            % starting today
            patchends = datenum([sv(1:3) 14 0 0]):1:t_interval(2);
        else % start is during the day locally
            % Patches begin at hour 2, starting w/ the next one (which may be today 
            % or tomorrow, depending on where t_interval(1) falls rel. midnight)
            % Patches end at hour 14, ditto, plus ditto comment above about end of last patch
            if shour <= 2 % it's after midnight, start and end patches begin today
                patchstarts = datenum([sv(1:3) 2 0 0]):1:t_interval(2);
                patchends = datenum([sv(1:3) 14 0 0]):1:t_interval(2);
            else % it's before midnight, start and end patches begin tomorrow
                patchstarts = datenum([sv(1:2) sv(3)+1 2 0 0]):1:t_interval(2);
                patchends = datenum([sv(1:2) sv(3)+1 14 0 0]):1:t_interval(2);
            end
        end
        if patchstarts(end) > patchends(end) % if there isn't an end for the last patchstart
            patchends = [patchends t_interval(2)];
        end

        % Create moving averages and NaNs here
        keys = {'P', 'T', 'S', 'DO'};
        vals = {P, T, S, O};
        Y_Map = Map(keys, vals);
        [twNaNs, PTSDOwNaNs, t_av, PTSDOav] = ...
            create_moving_av(t_interval, t, Y_Map, av_win, av_thresh);
        PwNaNs = PTSDOwNaNs('P');
        P_av = PTSDOav('P'); 
        TwNaNs = PTSDOwNaNs('T');
        T_av = PTSDOav('T'); 
        SwNaNs = PTSDOwNaNs('S');
        S_av = PTSDOav('S');
        OwNaNs = PTSDOwNaNs('DO');
        O_av = PTSDOav('DO');

        % Dissolved Oxygen time series (plotted first, despite being at the bottom of the "stack" of
        % time series plots, in order to cover the ticks otherwise displayed along this plot's top border)
        axes('position', [tserieslp tseriesbp(4) tseriesw tseriesh]);
        safe_t = t;
        if isnan(nanvar(O)) || nanvar(O) == 0
            t = t_interval(1):1/4/24:t_interval(2);
            O = floor(t/t(end));
            plotstr = 'plot(twNaNs, O, ''.b'', ''MarkerSize'', 1);';
        else
            plotstr = 'plot(twNaNs, OwNaNs, ''b'');';
        end
        addpatches(patchstarts, patchends, min(O), max(O));
        eval(plotstr);
        t = safe_t;
        % Note: plot has already been "held" by call to addpatches
        plot(t_av, O_av, '--c', 'LineWidth', 1.5);
        axis tight;
        xlabel('(UTC)');
        ylabel('Oxygen (mg/L)');
        Omin = min(O); Omax = max(O); Orange = Omax - Omin;
        Otickvals = Omin + Orange*[1 2 3 4 5]/6;
        Oticklabels = {sprintf('%4.1f', Otickvals(1)), '',...
                       sprintf('%4.1f', Otickvals(3)), '',...
                       sprintf('%4.1f', Otickvals(5))};
        xticks = ceil(t_interval(1)):1:t_interval(2);
        xticklabels = datestr(xticks, 'mm/dd');
        set(gca, 'xtickmode', 'manual',...
                 'TickDir', 'out',...
                 'TickLength', [0.01 0.01],...
                 'XLim', t_interval, 'XTick', xticks, 'XTickLabel', xticklabels,...
                 'ytick', Otickvals,...
                 'yticklabel', Oticklabels,...
                 'box', 'on');

        % Pressure time series
        axes('position', [tserieslp tseriesbp(1) tseriesw tseriesh]);
        safe_t = t;
        if isnan(nanvar(P)) || nanvar(P) == 0
            noP = true;
            t = t_interval(1):1/4/24:t_interval(2);
            P = floor(t/t(end));
            plotstr = 'plot(twNaNs, P, ''.k'', ''MarkerSize'', 1);';
        else
            plotstr = 'plot(twNaNs, PwNaNs, ''k'');';
            noP = false;
        end
        Pmin = min(P); Pmax = max(P);
        if strcmp(siteID(end), 'F') && ~noP
            Pmin = Pmin - 3; Pmax = Pmax + 3;
        end        
        addpatches(patchstarts, patchends, Pmin, Pmax);
        eval(plotstr);
        t = safe_t;
        % Note: plot has already been "held" by call to addpatches
        plot(t_av, P_av, '--k', 'LineWidth', 1.5);
        temp = text('Title');
        title(sprintf('%s %s (%s, %s)', temp{:}), 'FontSize', 14);
        ylabel('Pressure (dbar)');
        Prange = Pmax - Pmin;
        Ptickvals = Pmin + Prange*[1 2 3 4 5]/6;
        Pticklabels = {sprintf('%4.1f', Ptickvals(1)), '',...
                       sprintf('%4.1f', Ptickvals(3)), '',...
                       sprintf('%4.1f', Ptickvals(5))};
        set(gca, 'TickDir', 'out',...
                   'TickLength', [0.01 0.01],...
                   'xtick', [],...
                   'xlim', t_interval,...
                   'ylim', [Pmin Pmax],...
                   'ytick', Ptickvals,...
                   'yticklabel', Pticklabels,...
                   'box', 'on');
        axis tight;


        % Temperature time series
        axes('position', [tserieslp tseriesbp(2) tseriesw tseriesh]);
        safe_t = t;
        if isnan(nanvar(T)) || nanvar(T) == 0
            t = t_interval(1):1/4/24:t_interval(2);
            T = floor(t/t(end));
            plotstr = 'plot(twNaNs, T, ''.k'', ''MarkerSize'', 1);';
        else
            plotstr = 'plot(twNaNs, TwNaNs, ''r'');';        
        end
        addpatches(patchstarts, patchends, min(T), max(T));
        eval(plotstr);
        t = safe_t;
        % Note: plot has already been "held" by call to addpatches
        plot(t_av, T_av, '--m', 'LineWidth', 1.5);
        axis tight;
        ylabel('Temperature (\circC)');
        Tmin = min(T); Tmax = max(T); Trange = Tmax - Tmin;
        Ttickvals = Tmin + Trange*[1 2 3 4 5]/6;
        Tticklabels = {sprintf('%4.1f', Ttickvals(1)), '',...
                       sprintf('%4.1f', Ttickvals(3)), '',...
                       sprintf('%4.1f', Ttickvals(5))};
        set(gca, 'TickDir', 'out',...
                   'TickLength', [0.01 0.01],...
                   'xtick', [],...
                   'xlim', t_interval,...
                   'ytick', Ttickvals,...
                   'yticklabel', Tticklabels,...
                   'box', 'on')

        % Salinity time series
        axes('position', [tserieslp tseriesbp(3) tseriesw tseriesh]);
        safe_t = t;
        if isnan(nanvar(S)) || nanvar(S) == 0
            t = t_interval(1):1/4/24:t_interval(2);
            S = floor(t/t(end));
            plotstr = 'h = plot(twNaNs, S, ''.'', ''MarkerSize'', 1);';
            skip_contour = true;
        else
            plotstr = 'h = plot(twNaNs, SwNaNs);';        
        end
        addpatches(patchstarts, patchends, min(S), max(S))
        eval(plotstr);
        set(h, 'color', [0.87 0.49 0]);
        t = safe_t;
        % Note: plot has already been "held" by call to addpatches
        h = plot(t_av, S_av, '--', 'LineWidth', 1.5);
        set(h, 'color', [0.9 0.6 0]);
        axis tight;
        ylabel('Salinity (PSU)');
        Smin = min(S); Smax = max(S); Srange = Smax - Smin;
        Stickvals = Smin + Srange*[1 2 3 4 5]/6;
        Sticklabels = {sprintf('%4.1f', Stickvals(1)), '',...
                       sprintf('%4.1f', Stickvals(3)), '',...
                       sprintf('%4.1f', Stickvals(5))};
        set(gca, 'TickDir', 'out',...
                   'TickLength', [0.01 0.01],...
                   'xtick', [],...
                   'xlim', t_interval,...
                   'ytick', Stickvals,...
                   'yticklabel', Sticklabels,...
                   'box', 'on');

        % (Water) Density(S,T) and CDF(S,T) contours
        axes('position', [tserieslp 0.115 0.3 0.175], 'color', 'k', 'layer', 'bottom');
        hold;
        nS = 20; nT = 20;
        levels = [0.05 0.25 0.5 0.75 0.95];
        SM = max(S); Sm = min(S); Srange = SM - Sm;
        TM = max(T); Tm = min(T); Trange = TM - Tm;
        tol = 0.2;
        if ~skip_contour
            rt_2dhist(gca, S, T,...
                      [(Sm - tol*Srange) (SM + tol*Srange) nS],...
                      [(Tm - tol*Trange) (TM + tol*Trange) nT],...
                      levels);
            Dofdata = sw_dens0(S, T) - 1000; 
            range = max(Dofdata)-min(Dofdata);
            incrs = [0.25 0.5 1.0 1.5]; 
            ranges = 5*incrs;
            dr = abs(ranges-range);
            incr = incrs(dr == min(dr));
            Sbins = linspace(10, 35, nS);
            Tbins = linspace(0, 25, nT);
            [SxS, TxT] = meshgrid(Sbins, Tbins);
            D = sw_dens0(SxS, TxT) - 1000;
            minD = round2(min(min(Dofdata)), incr, 'floor'); 
            maxD = round2(max(max(Dofdata)), incr, 'ceil');
            [C, hc] = contour(SxS, TxT, D, minD:incr:maxD);
            set(gca, 'TickDir', 'out',...
                     'TickLength', [0.01 0.01]);
            set(hc, 'LineColor', 'w');
            clabelH = clabel(C, hc, 'LabelSpacing', 144);
            set(clabelH(:), 'FontSize', 8,...
                            'Color', 'w');
            set(gca, 'Box', 'on');
            title('Density (kg/m^3)');
            xlabel('Salinity (PSU)');
            ylabel('Temperature (\circC)');
            tmp = get(gca, {'XLim', 'YLim', 'XTick', 'YTick', 'XTickLabel', 'YTickLabel'});
            [xlims, ylims, xticks, yticks, xtlabs, ytlabs] = tmp{:};
            cmap = colormap; 
            cmapres = length(cmap(:,1)); 
            cs = colorscale([1 cmapres], [levels(1) levels(end)], 0.05, 'vertical',...
                 'position', [0.49 0.115 0.01 0.175], 'tickdir', 'out');
            axis tight;
            csticklabels = cellstr(num2str(levels(:)));
            csticks = cmapres*levels; 
            set(cs, 'YTick', csticks, 'YTickLabel', csticklabels, 'FontSize', 8, 'Box', 'off');
            ylabelH = get(cs, 'YLabel');
            pos = get(ylabelH, 'Position');
            set(ylabelH, 'String', 'Probability', 'FontSize', 9, 'Position', [-3 pos(2:3)]);
        end

        % Dissolved Oxygen(S,T) scatter plot
        site = text('Title');
%        mch = strcmp(site{1}, 'MCH01BR');
        if ~isempty(O) && any(~isnan(O))% || mch
            axes('position', [0.6 0.115 0.3 0.175]);
            if ~skip_contour
                H = scatter(S, T, 'CData', -O, 'sizedata', 15);
                hold;
                caxis([-7 -4]);
                set(gca, 'TickDir', 'out', 'Box', 'on', 'YaxisLocation', 'right',...
                       'XLim', xlims, 'XTick', xticks, 'XTickLabel', xtlabs,...
                       'YLim', ylims, 'YTick', yticks, 'YTickLabel', ytlabs);
                [C, h] = contour(SxS, TxT, D, minD:incr:maxD);
                set(h, 'LineColor', [0.35 0.35 0.35]);
                clabelH = clabel(C, h, 'LabelSpacing', 144);
                set(clabelH(:), 'FontSize', 8,...
                                'Color', [0.35, 0.35, 0.35]);
                title('Oxygen (DO) (mg/L)');
                xlabel('Salinity (PSU)');
                ylabel('Temperature (\circC)');
                levels = [4 4.75 5.5 6.25 7];
                levellabels = {'4' '4.75' '5.5' '6.25' '7'};
                yincr = (levels(1)-levels(end))/(cmapres-1);        
                cs = colorscale([cmapres 1], [levels(1) levels(end)], abs(yincr), 'vertical',... 
                                'position', [0.5 0.115 0.01 0.175], 'YDir', 'reverse', 'tickdir', 'out',...
                                'YaxisLocation', 'right', 'YTick', 0.5:cmapres/4:cmapres+0.5,...
                                'YTickLabel', levellabels, 'FontSize', 8, 'Box', 'off');
                ylabelH = get(cs, 'YLabel');
                pos = get(ylabelH, 'Position');
                set(ylabelH, 'String', 'DO (mg/L)', 'FontSize', 9, 'Position', [5 pos(2:3)]);
            end
        end

        % Caption
        [caption, r] = strtok(text('Caption'), ';');
        caption = strrep(caption, '__', num2str(24*av_win));
        [c, r] = strtok(r, ';');
        caption = [caption char(10) c];
        if ~isempty(O) && any(~isnan(O))
            [c, r] = strtok(r, ';');
            caption = [caption char(10) c];
        end
        annotation('textbox', [0.05 0.025 0.9 0.05], 'String', caption,...
                   'LineStyle', 'none', 'fontsize', 8);
    catch PE
        figure_handle = PE;
    end
end

function addpatches(ps, pe, ymin, ymax)
    % Function to draw the nighttime-indicating patches
    pcol = [0.85 0.85 0.85];
    for i=1:length(ps)
        patch([ps(i) pe(i) pe(i) ps(i)],...
              [ymin ymin ymax ymax],...
              pcol, 'edgecolor', 'none');
    end
    hold;
end
        
function [twNaNs, YwNaNs, t_av, Y_av] = create_moving_av(t_interval, t, Y, window, threshold)
    % Function to create moving averages, as well as NaN missing data
    import containers.*

    Ykeys = keys(Y);
    Yvals = values(Y);
    dt = 1/24/4; % 15 minutes in units of days
    % Create NaN filled arrays
    twNaNs = t_interval(1):dt:t_interval(2);
    YwNaNstmp = nans([Y.Count length(twNaNs)]); % Array storage for YwNaNs
    if ~isempty(Y('T'))
        for i=1:length(twNaNs)
            tmp = find((t-dt/2 <= twNaNs(i)) & (twNaNs(i) <= t+dt/2));
            if ~isempty(tmp)
                twNaNs(i) = t(tmp(1));
                for k=1:int8(Y.Count)
                    Y_tmp = Yvals{k};
                    if ~isempty(Y_tmp)
                        YwNaNstmp(k,i) = Y_tmp(tmp(1));
                    end
                end
            end
        end
    end
    % Convert YwNaNstmp to Map object
    YwNaNs = Map();
    for k=1:int8(Y.Count)
        YwNaNs(Ykeys{k}) = YwNaNstmp(k,:);
    end
    
    % Create moving averages
    half_win = window/2;
    t_av = (t_interval(1)+half_win):dt:(t_interval(2)-half_win);
    Y_avtmp = nans([Y.Count length(t_av)]); % Array storage for Y_av
    if ~isempty(Y('T'))
        for i=1:length(t_av)
            tmp = find((t_av(i)-half_win <= t) & (t <= t_av(i)+half_win));
            if length(tmp) >= threshold
                for k=1:int8(Y.Count)
                    Y_tmp = Yvals{k};
                    if ~isempty(Y_tmp)
                        Y_avtmp(k,i) = nanmean(Y_tmp(tmp));
                    end
                end
            end
        end
    end
    % Convert Y_avtmp to Map object
    Y_av = Map();
    Ykeys = keys(Y);
    for k=1:int8(Y.Count)
        Y_av(Ykeys{k}) = Y_avtmp(k,:);
    end
end